options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(1,4))
      drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
      drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
      drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
      drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}

normal distribution

ex1.stan

normal distribution

data{
  int N;
  vector[N] y;
}

parameters{
  real m;
  real<lower=0> s;
}

model{
  y~normal(m,s);
}
mdl=cmdstan_model('./ex1.stan')
y=rnorm(10,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.52    1.06    1.74    1.86    2.55    4.14
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')  

data=list(N=length(y),y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -13.6463 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       10      -5.59216   7.78142e-05   3.67445e-05           1           1       13    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.3 seconds.
mle
##  variable estimate
##      lp__    -5.59
##      m        1.86
##      s        1.06
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
mcmc
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -6.57  -6.26 1.03 0.76 -8.58 -5.57 1.00      715      752
##      m     1.87   1.87 0.43 0.41  1.17  2.56 1.01      697      714
##      s     1.30   1.24 0.36 0.28  0.86  1.93 1.00     1295      960
mcmc$metadata()$stan_variables
## [1] "lp__" "m"    "s"
mcmc$metadata()$model_params
## [1] "lp__" "m"    "s"
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -6.57  -6.26 1.03 0.76 -8.58 -5.57 1.00      715      752
##      m     1.87   1.87 0.43 0.41  1.17  2.56 1.01      697      714
##      s     1.30   1.24 0.36 0.28  0.86  1.93 1.00     1295      960

poisson distribution

y=rpois(20,3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    2.00    3.00    3.05    4.00    7.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=length(y),y=y)

ex2-1.stan

in case fit poisson dist. to normal dist.

data{
  int N;
  vector[N] y;
}

parameters{
  real<lower=0> m;
  real<lower=0> s;
}

model{
  y~normal(m,s);
}

generated quantities{
  vector[N] y1;
  for(i in 1:N)
    y1[i]=normal_rng(m,s);
}
# when fitting poisson dist. sample to normal dist.
mdl=cmdstan_model('./ex2-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -34.2029 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       11       -17.643    0.00347767    0.00136638      0.9986      0.9986       24    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -17.64
##    m          3.05
##    s          1.47
##    y1[1]      1.29
##    y1[2]      1.93
##    y1[3]      3.15
##    y1[4]      2.82
##    y1[5]      3.54
##    y1[6]      3.69
##    y1[7]      2.17
##    y1[8]      3.85
##    y1[9]      2.87
##    y1[10]     2.60
##    y1[11]     2.47
##    y1[12]     0.96
##    y1[13]     5.63
##    y1[14]     3.03
##    y1[15]     1.84
##    y1[16]     3.32
##    y1[17]     3.20
##    y1[18]     0.70
##    y1[19]     3.26
##    y1[20]     2.78
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -17.28 -16.90 1.21 0.81 -19.60 -16.19 1.00      852      949
##    m        3.05   3.04 0.38 0.36   2.42   3.67 1.00     1747     1065
##    s        1.63   1.59 0.31 0.28   1.21   2.16 1.01      825      802
##    y1[1]    3.12   3.10 1.74 1.61   0.37   6.04 1.00     1964     1834
##    y1[2]    3.06   3.07 1.72 1.63   0.24   5.96 1.00     1754     1759
##    y1[3]    3.08   3.11 1.69 1.55   0.34   5.82 1.00     2093     1875
##    y1[4]    3.00   3.06 1.70 1.68   0.16   5.70 1.00     1785     1801
##    y1[5]    2.98   2.98 1.69 1.61   0.16   5.68 1.00     2021     1870
##    y1[6]    3.09   3.08 1.70 1.63   0.32   5.90 1.00     2071     2027
##    y1[7]    3.03   2.95 1.71 1.63   0.34   5.92 1.00     1993     2054
##    y1[8]    3.01   3.05 1.70 1.59   0.16   5.86 1.00     1532     1593
##    y1[9]    2.99   2.96 1.69 1.61   0.25   5.88 1.00     1952     1931
##    y1[10]   3.07   3.10 1.74 1.69   0.23   5.95 1.00     1594     1927
##    y1[11]   3.03   3.03 1.72 1.68   0.22   5.79 1.00     2015     1995
##    y1[12]   3.08   3.08 1.75 1.73   0.20   5.90 1.00     2074     1890
##    y1[13]   3.02   3.05 1.65 1.56   0.28   5.70 1.00     2050     1970
##    y1[14]   2.98   2.98 1.70 1.69   0.18   5.78 1.00     2034     1895
##    y1[15]   3.08   3.08 1.65 1.59   0.43   5.73 1.00     1897     1916
##    y1[16]   3.02   2.98 1.73 1.67   0.26   5.80 1.00     1987     1840
##    y1[17]   2.99   3.03 1.66 1.56   0.24   5.66 1.00     2021     1849
##    y1[18]   3.07   3.07 1.70 1.67   0.30   5.88 1.00     1857     1793
##    y1[19]   3.07   3.05 1.71 1.70   0.29   5.86 1.00     2049     2057
##    y1[20]   3.09   3.09 1.72 1.56   0.22   5.91 1.00     1776     1725

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1')     

ex2-2.stan

fit to poisson dist.

data{
  int N;
  array[N] int y;
}

parameters{
  real<lower=0> l;
}

model{
  y~poisson(l);
}

generated quantities{
  array[N] int y1;
  for(i in 1:N)
    y1[i]=poisson_rng(l);
}
mdl=cmdstan_model('./ex2-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -16.1399 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6       7.02364   0.000192372   3.67643e-05           1           1        7    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__       7.02
##    l          3.05
##    y1[1]      4.00
##    y1[2]      2.00
##    y1[3]      4.00
##    y1[4]      3.00
##    y1[5]      2.00
##    y1[6]      1.00
##    y1[7]      3.00
##    y1[8]      1.00
##    y1[9]      1.00
##    y1[10]     4.00
##    y1[11]     7.00
##    y1[12]     2.00
##    y1[13]     1.00
##    y1[14]     1.00
##    y1[15]     2.00
##    y1[16]     4.00
##    y1[17]     4.00
##    y1[18]     3.00
##    y1[19]     4.00
##    y1[20]     0.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable mean median   sd  mad   q5  q95 rhat ess_bulk ess_tail
##    lp__   7.66   7.93 0.69 0.30 6.29 8.15 1.00     1032     1449
##    l      3.10   3.10 0.39 0.38 2.48 3.73 1.00      817     1208
##    y1[1]  3.09   3.00 1.81 1.48 1.00 6.00 1.00     1805     1596
##    y1[2]  3.11   3.00 1.82 1.48 0.95 6.00 1.00     1666     1870
##    y1[3]  3.11   3.00 1.75 1.48 1.00 6.00 1.00     1922     1633
##    y1[4]  3.03   3.00 1.80 1.48 0.00 6.00 1.00     1877     1939
##    y1[5]  3.10   3.00 1.81 1.48 0.00 6.00 1.00     1998     1622
##    y1[6]  3.11   3.00 1.83 1.48 0.00 6.00 1.00     1932     1788
##    y1[7]  3.06   3.00 1.79 1.48 1.00 6.00 1.00     1786     1868
##    y1[8]  3.01   3.00 1.78 1.48 0.00 6.00 1.00     1835     1744
##    y1[9]  3.12   3.00 1.76 1.48 1.00 6.00 1.00     1917     1960
##    y1[10] 3.12   3.00 1.77 1.48 0.00 6.00 1.00     1741     1947
##    y1[11] 3.04   3.00 1.75 1.48 0.00 6.00 1.00     1973     1988
##    y1[12] 3.04   3.00 1.77 1.48 0.00 6.00 1.00     1951     1832
##    y1[13] 3.10   3.00 1.78 1.48 1.00 6.00 1.00     1987     2020
##    y1[14] 3.01   3.00 1.80 1.48 0.00 6.00 1.00     2104     1909
##    y1[15] 3.14   3.00 1.82 1.48 1.00 6.00 1.00     1889     1495
##    y1[16] 3.08   3.00 1.85 1.48 1.00 7.00 1.00     1934     1799
##    y1[17] 3.10   3.00 1.79 1.48 1.00 6.00 1.00     1928     1994
##    y1[18] 3.08   3.00 1.78 1.48 1.00 6.00 1.00     2085     1995
##    y1[19] 3.14   3.00 1.79 1.48 1.00 6.00 1.00     1730     1877
##    y1[20] 3.12   3.00 1.81 1.48 0.00 6.00 1.00     1656     1811

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

difference test

ex3.stan

mean difference of 2 normal distributions

data{
  int N;
  vector[N] a;
  vector[N] b;
}
parameters{
  real ma;
  real<lower=0> sa;
  real mb;
  real<lower=0> sb;
}
model{
  a~normal(ma,sa);
  b~normal(mb,sb);
}
generated quantities{
  real d;
  d=mb-ma;
}
n=20
a=rnorm(n,10,1)
b=rnorm(n,11,2)
data=list(N=n,a=a,b=b)

mdl=cmdstan_model('./ex3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -29116.3 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       19      -32.5812   0.000187185   0.000717354           1           1       52    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -32.58
##      ma      10.14
##      sa       1.08
##      mb      10.71
##      sb       1.74
##      d        0.56
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -34.01 -33.66 1.49 1.25 -36.93 -32.32 1.00      772     1202
##      ma    10.13  10.13 0.26 0.26   9.72  10.54 1.00     1290     1090
##      sa     1.19   1.16 0.22 0.20   0.89   1.58 1.00     2567     1394
##      mb    10.73  10.73 0.43 0.41  10.02  11.44 1.01      824      769
##      sb     1.91   1.86 0.34 0.31   1.45   2.56 1.01     1953     1444
##      d      0.59   0.59 0.50 0.50  -0.26   1.42 1.01      899      938

d=mcmc$draws('d')
d
## # A draws_array: 500 iterations, 4 chains, and 1 variables
## , , variable = d
## 
##          chain
## iteration    1     2    3    4
##         1 0.18 -0.23 0.96 0.22
##         2 0.30  1.16 0.96 0.58
##         3 0.62  1.09 1.10 0.53
##         4 0.50  0.97 0.91 0.66
##         5 0.52  1.08 0.95 1.13
## 
## # ... with 495 more iterations
mean(d>0)
## [1] 0.883

normal regression

ex4-1.stan

single normal regression

data{
  int N;
  vector[N] x;
  vector[N] y;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}
# make prediction with R
n=20
x=runif(n,0,100)
y=rnorm(n,x*3+10,10)
par(mfrow=c(1,1))
plot(x,y)

data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex4-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -824269 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpSt0wv7/model-a38248be373.stan', line 14, column 2 to column 22) 
## Error evaluating model log probability: Non-finite gradient. 
##       42      -53.5663    0.00362798   0.000931153      0.9809      0.9809      108    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -53.57
##      b0      -2.64
##      b1       3.20
##      s        8.83
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.2 seconds.
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -53.06 -52.69 1.42 1.12 -55.83 -51.55 1.01      319      303
##      b0    -2.25  -2.61 6.21 5.68 -12.11   8.26 1.02      175      104
##      b1     3.20   3.20 0.09 0.08   3.05   3.33 1.02      198      149
##      s     10.08   9.81 1.83 1.63   7.55  13.46 1.00      814      917

b0=mcmc$draws('b0')
b1=mcmc$draws('b1')
b=tibble(b0=c(b0),b1=c(b1)) |> sample_n(100)
x1=runif(nrow(b),min(x),max(x))
xy=tibble(x=x1,y=b$b0+b$b1*x1)
par(mfrow=c(1,1))
plot(xy)

ex4-2.stan

single normal regression, prediction from new explanatory variable

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}

generated quantities{
  vector[N1] m;
  vector[N1] y1;
  for(i in 1:N1){
    m[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m[i],s);
  }
}
# make prediction with stan
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex4-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -18297.4 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       33      -53.5663     0.0672733    0.00109384           1           1       75    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -53.57
##    b0        -2.64
##    b1         3.20
##    s          8.83
##    m[1]      -1.79
##    m[2]      33.55
##    m[3]      68.89
##    m[4]     104.23
##    m[5]     139.56
##    m[6]     174.90
##    m[7]     210.24
##    m[8]     245.58
##    m[9]     280.92
##    m[10]    316.26
##    y1[1]     -1.50
##    y1[2]     33.50
##    y1[3]     71.09
##    y1[4]    108.73
##    y1[5]    160.00
##    y1[6]    168.09
##    y1[7]    202.29
##    y1[8]    248.71
##    y1[9]    272.47
##    y1[10]   325.25
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 0.3 seconds.
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "s"    "m"    "y1"
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median    sd   mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -53.04 -52.62  1.38  1.10 -55.70 -51.55 1.01      330      423
##    b0      -1.62  -1.96  5.94  5.56 -11.02   9.34 1.02      199      164
##    b1       3.19   3.20  0.08  0.08   3.05   3.33 1.02      226      273
##    s        9.98   9.73  1.81  1.71   7.53  13.34 1.00      864      787
##    m[1]    -0.78  -1.12  5.92  5.54 -10.14  10.15 1.02      199      164
##    m[2]    34.42  34.14  5.08  4.67  26.41  43.77 1.02      201      125
##    m[3]    69.62  69.42  4.28  3.95  62.88  77.37 1.03      206      177
##    m[4]   104.82 104.74  3.54  3.31  99.32 111.07 1.02      218      226
##    m[5]   140.02 139.98  2.91  2.77 135.50 145.01 1.02      256      363
##    m[6]   175.22 175.17  2.46  2.31 171.31 179.37 1.00      396      532
##    m[7]   210.43 210.42  2.31  2.20 206.67 214.39 1.01     1106     1062
##    m[8]   245.63 245.66  2.51  2.42 241.53 249.81 1.01     2122     1123
##    m[9]   280.83 280.87  2.99  2.87 275.75 285.63 1.01     1254     1084
##    m[10]  316.03 316.11  3.65  3.63 309.75 321.87 1.01      688     1012
##    y1[1]   -0.73  -0.64 11.66 11.82 -19.35  19.11 1.00      613     1204
##    y1[2]   34.33  33.90 11.24 10.84  15.96  53.16 1.00      670     1026
##    y1[3]   69.82  69.86 11.16 10.59  50.96  87.72 1.00      978     1374
##    y1[4]  105.03 105.14 10.77  9.97  87.33 122.45 1.00      865     1380
##    y1[5]  139.98 140.20 10.75 10.60 122.25 157.46 1.00     1722     1549
##    y1[6]  175.51 175.70 10.56 10.32 158.27 192.08 1.00     1797     1835
##    y1[7]  210.27 210.08 10.40 10.20 193.35 227.49 1.00     1970     1869
##    y1[8]  245.72 245.83 10.12 10.16 229.02 261.62 1.00     1910     1870
##    y1[9]  281.00 280.89 10.58  9.76 262.83 298.58 1.00     1798     1646
##    y1[10] 315.80 315.58 10.91 10.72 297.93 334.03 1.00     1767     1525

m=mcmc$draws('m')
summary(m)
## # A tibble: 10 × 10
##    variable    mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
##    <chr>      <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
##  1 m[1]      -0.779  -1.12  5.92  5.54 -10.1  10.2  1.02     200.     164.
##  2 m[2]      34.4    34.1   5.08  4.67  26.4  43.8  1.02     202.     126.
##  3 m[3]      69.6    69.4   4.28  3.95  62.9  77.4  1.03     206.     177.
##  4 m[4]     105.    105.    3.54  3.31  99.3 111.   1.02     218.     226.
##  5 m[5]     140.    140.    2.91  2.77 135.  145.   1.02     257.     364.
##  6 m[6]     175.    175.    2.46  2.31 171.  179.   1.00     396.     533.
##  7 m[7]     210.    210.    2.31  2.20 207.  214.   1.01    1106.    1062.
##  8 m[8]     246.    246.    2.51  2.42 242.  250.   1.01    2122.    1123.
##  9 m[9]     281.    281.    2.99  2.87 276.  286.   1.01    1255.    1085.
## 10 m[10]    316.    316.    3.65  3.63 310.  322.   1.01     688.    1013.
smm=summary(m)

y1=mcmc$draws('y1')
summary(y1)
## # A tibble: 10 × 10
##    variable    mean  median    sd   mad    q5   q95  rhat ess_bulk ess_tail
##    <chr>      <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
##  1 y1[1]     -0.730  -0.645  11.7 11.8  -19.4  19.1 1.00      614.    1204.
##  2 y1[2]     34.3    33.9    11.2 10.8   16.0  53.2 1.00      670.    1026.
##  3 y1[3]     69.8    69.9    11.2 10.6   51.0  87.7 1.00      979.    1375.
##  4 y1[4]    105.    105.     10.8  9.97  87.3 122.  1.00      865.    1381.
##  5 y1[5]    140.    140.     10.8 10.6  122.  157.  1.00     1722.    1550.
##  6 y1[6]    176.    176.     10.6 10.3  158.  192.  0.999    1797.    1836.
##  7 y1[7]    210.    210.     10.4 10.2  193.  227.  1.00     1971.    1870.
##  8 y1[8]    246.    246.     10.1 10.2  229.  262.  1.00     1911.    1871.
##  9 y1[9]    281.    281.     10.6  9.76 263.  299.  1.00     1798.    1647.
## 10 y1[10]   316.    316.     10.9 10.7  298.  334.  1.00     1768.    1526.
smy=summary(y1)

xy=tibble(x=x1,m=smm$median,ml5=smm$q5,mu5=smm$q95,yl5=smy$q5,yu5=smy$q95)

par(mfrow=c(1,1))
xlim=c(min(x),max(x))
ylim=c(min(y),max(y))
plot(x,y,
     xlim=xlim,ylim=ylim, xlab='x',ylab='y',col='red')
par(new=T)
plot(xy$x,xy$m,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='black')
par(new=T)
plot(xy$x,xy$ml5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$mu5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$yl5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
par(new=T)
plot(xy$x,xy$yu5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=m),xy,col='black')+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')

ex4-3.stan

single normal regression, prediction from original explanatory variable

data{
  int N;
  vector[N] x;
  vector[N] y;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}

generated quantities{
  vector[N] m;
  vector[N] y1;
  for(i in 1:N){
    m[i]=b0+b1*x[i];
    y1[i]=normal_rng(m[i],s);
  }
}
data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex4-3.stan')

mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median    sd   mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -53.04 -52.68  1.33  1.10 -55.64 -51.58 1.01      364      350
##    b0      -2.95  -2.95  6.14  5.92 -12.39   7.56 1.03      151      186
##    b1       3.21   3.21  0.09  0.08   3.07   3.35 1.03      174      315
##    s       10.11   9.92  1.76  1.67   7.62  13.49 1.01      781      690
##    m[1]   188.92 188.86  2.35  2.26 185.04 192.76 1.01      608      688
##    m[2]   309.72 309.62  3.63  3.54 304.00 315.64 1.01      619     1248
##    m[3]   175.51 175.48  2.45  2.27 171.50 179.59 1.01      382      650
##    m[4]   262.90 262.89  2.76  2.70 258.47 267.47 1.01     1720     1664
##    m[5]   176.45 176.41  2.44  2.26 172.46 180.52 1.01      393      649
##    m[6]    71.79  71.74  4.32  4.01  65.10  79.27 1.03      153      204
##    m[7]   121.32 121.24  3.26  2.98 116.25 126.91 1.03      171      278
##    m[8]   110.26 110.20  3.48  3.20 104.84 116.27 1.03      163      251
##    m[9]   111.67 111.61  3.45  3.17 106.29 117.62 1.03      163      248
##    m[10]  239.56 239.51  2.47  2.49 235.59 243.72 1.01     2227     1609
##    m[11]  194.98 194.93  2.32  2.25 191.16 198.77 1.01      722      696
##    m[12]  306.21 306.10  3.55  3.46 300.61 312.01 1.01      651     1244
##    m[13]  309.30 309.20  3.62  3.53 303.59 315.20 1.01      623     1248
##    m[14]  219.51 219.51  2.33  2.34 215.68 223.28 1.00     1694     1571
##    m[15]   -2.10  -2.10  6.12  5.88 -11.49   8.38 1.03      151      186
##    m[16]  153.46 153.45  2.71  2.50 149.05 158.06 1.02      232      426
##    m[17]  308.54 308.45  3.60  3.51 302.85 314.43 1.01      630     1248
##    m[18]  316.48 316.35  3.77  3.68 310.49 322.70 1.01      569     1102
##    m[19]  265.53 265.51  2.80  2.75 261.03 270.14 1.01     1650     1690
##    m[20]  179.65 179.59  2.41  2.24 175.71 183.68 1.01      438      665
##    y1[1]  188.67 188.62 10.20  9.84 172.14 205.54 1.00     1622     1691
##    y1[2]  309.93 310.02 10.94 10.22 290.90 328.43 1.00     1737     1957
##    y1[3]  175.71 175.81 10.27  9.70 158.54 192.63 1.00     1869     1870
##    y1[4]  262.53 262.67 10.49 10.21 244.88 279.56 1.00     1993     1880
##    y1[5]  176.65 176.59 10.80 10.64 159.09 193.76 1.00     1902     1972
##    y1[6]   71.86  71.81 11.19 10.92  53.75  90.43 1.00     1008      903
##    y1[7]  121.49 121.32 10.83 10.27 103.52 139.28 1.00     1499     1463
##    y1[8]  110.72 110.81 10.70 10.34  93.04 128.44 1.00     1164     1599
##    y1[9]  111.97 112.07 10.80 10.25  94.21 129.81 1.00     1240     1733
##    y1[10] 239.66 239.61 10.69 10.49 222.17 257.29 1.00     2093     2002
##    y1[11] 195.09 195.16 10.35 10.06 178.41 212.12 1.00     1949     1723
##    y1[12] 306.11 306.36 10.61 10.73 288.36 323.52 1.00     1884     1928
##    y1[13] 309.27 309.20 10.94 10.49 291.03 327.03 1.00     1539     1958
##    y1[14] 219.16 219.28 10.54  9.79 201.85 236.33 1.00     1791     1968
##    y1[15]  -2.49  -2.31 12.05 11.36 -22.20  17.87 1.01      647     1304
##    y1[16] 153.57 153.71 10.65 10.39 136.56 170.59 1.00     1736     1773
##    y1[17] 308.13 308.30 10.78 10.34 290.42 325.66 1.00     1770     1999
##    y1[18] 316.78 316.78 10.85 10.75 298.57 333.96 1.00     1432     1874
##    y1[19] 265.50 265.31 10.57 10.42 248.04 282.84 1.00     1980     1886
##    y1[20] 179.84 179.95 10.58 10.43 162.19 196.68 1.00     1675     1883

y1=mcmc$draws('y1')
smy=summary(y1)

par(mfrow=c(1,1))
plot(y,smy$median,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(y,smy$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(y-smy$median,xlab='obs.-prd.',main='residual')
density(y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

estimate correlation

ex4-41.stan

use covariance matrix and multinormal distribution

data{
  int N;
  array[N] vector[2] y;
}

parameters{
  vector[2] m;
  vector<lower=0>[2] s;
  real<lower=-1,upper=1> r;
}

transformed parameters{
  cov_matrix[2] cv;
  cv[1,1]=s[1]^2;
  cv[2,2]=s[2]^2;
  cv[1,2]=r*s[1]*s[2];
  cv[2,1]=r*s[1]*s[2];
}

model{
  y~multi_normal(m,cv);
}

ex4-42.stan

use covariance matrix and wishart distribution

data{
  int N;
  matrix[2,2] dp;
}

parameters{
  vector<lower=0>[2] s;
  real<lower=-1,upper=1> r;
}

transformed parameters{
  cov_matrix[2] cv;
  cv[1,1]=s[1]^2;
  cv[2,2]=s[2]^2;
  cv[1,2]=r*s[1]*s[2];
  cv[2,1]=r*s[1]*s[2];
}

model{
  dp~wishart(N-1,cv);
}

ex4-43.stan

use correlation matrix and wishart distribution

data{
  int N;
  int M;
  matrix[M,M] dp;
}

parameters{
  vector<lower=0>[M] s;
  corr_matrix[M] r;
}

transformed parameters{
  cov_matrix[M] cv;
  cv=quad_form_diag(r,s);
}

model{
  dp~wishart(N-1,cv);
}

ex4-44.stan use covariance matrix and multinormal distribution

data{
  int N;
  int K;
  array[N] vector[K] y;
}
parameters{
  vector[K] m;
  cov_matrix[K] cov;
}
model{
  y~multi_normal(m,cov);
}
library(MASS)
n=20
y=mvrnorm(n,c(0,0),matrix(c(1,0.7,0.7,1),2),2)
cov(y)
##       [,1]  [,2]
## [1,] 0.997 0.590
## [2,] 0.590 0.843
cor(y)
##       [,1]  [,2]
## [1,] 1.000 0.643
## [2,] 0.643 1.000
plot(y)

#estimate covariance
mdl=cmdstan_model('./ex4-41.stan')

data=list(N=n,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -146.414 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       19      -11.9003   3.45688e-05   3.28354e-05           1           1       22    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__      -11.90
##   m[1]       -0.09
##   m[2]       -0.01
##   s[1]        0.97
##   s[2]        0.89
##   r           0.64
##   cv[1,1]     0.95
##   cv[2,1]     0.56
##   cv[1,2]     0.56
##   cv[2,2]     0.80
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -15.86 -15.50 1.63 1.50 -19.04 -13.84 1.01      539     1060
##   m[1]     -0.09  -0.09 0.24 0.23  -0.49   0.31 1.00     1413     1004
##   m[2]     -0.02  -0.02 0.22 0.22  -0.37   0.34 1.00     1450     1105
##   s[1]      1.08   1.05 0.19 0.18   0.81   1.44 1.00     1638     1471
##   s[2]      0.99   0.97 0.17 0.15   0.74   1.31 1.00     1655     1570
##   r         0.59   0.61 0.15 0.14   0.29   0.80 1.00      592      761
##   cv[1,1]   1.19   1.11 0.43 0.38   0.66   2.07 1.00     1638     1471
##   cv[2,1]   0.65   0.59 0.31 0.26   0.24   1.24 1.00      718      895
##   cv[1,2]   0.65   0.59 0.31 0.26   0.24   1.24 1.00      718      895
##   cv[2,2]   1.00   0.93 0.36 0.29   0.55   1.71 1.00     1655     1570

#estimate correlation,covariance
#deviation product sum matirix folows Wishart dist.
#deviation product sum = covariance * n-1 or n

mdl=cmdstan_model('./ex4-42.stan')

data=list(N=n,dp=cov(y)*(n-1))

mle=mdl$optimize(data=data)
## Initial log joint probability = -691.138 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       15      -12.2798   0.000110696   0.000833983           1           1       18    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__      -12.28
##   s[1]        1.00
##   s[2]        0.92
##   r           0.64
##   cv[1,1]     1.00
##   cv[2,1]     0.59
##   cv[1,2]     0.59
##   cv[2,2]     0.84
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -15.05 -14.71 1.22 0.95 -17.51 -13.74 1.00      920     1093
##   s[1]      1.08   1.06 0.19 0.17   0.82   1.45 1.00     1351     1384
##   s[2]      0.99   0.97 0.17 0.16   0.75   1.30 1.00     1228     1350
##   r         0.60   0.62 0.15 0.14   0.34   0.80 1.01      469      713
##   cv[1,1]   1.21   1.12 0.46 0.36   0.68   2.09 1.00     1351     1384
##   cv[2,1]   0.67   0.61 0.32 0.28   0.28   1.27 1.01      536      930
##   cv[1,2]   0.67   0.61 0.32 0.28   0.28   1.27 1.01      536      930
##   cv[2,2]   1.01   0.94 0.37 0.30   0.56   1.70 1.00     1229     1350

#estimate correlation,covariance
mdl=cmdstan_model('./ex4-43.stan')

m=2
data=list(N=n,M=m,dp=cov(y)*(n-1))

mle=mdl$optimize(data=data)
## Initial log joint probability = -50.0993 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       14      -12.2798   0.000249903   0.000125081           1           1       24    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__      -12.28
##   s[1]        1.00
##   s[2]        0.92
##   r[1,1]      1.00
##   r[2,1]      0.64
##   r[1,2]      0.64
##   r[2,2]      1.00
##   cv[1,1]     1.00
##   cv[2,1]     0.59
##   cv[1,2]     0.59
##   cv[2,2]     0.84
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -14.37 -14.02 1.27 0.97 -16.86 -13.04 1.00      796      786
##   s[1]      1.08   1.05 0.18 0.17   0.82   1.41 1.00     1360     1029
##   s[2]      0.99   0.97 0.17 0.16   0.76   1.32 1.00     1355     1211
##   r[1,1]    1.00   1.00 0.00 0.00   1.00   1.00   NA       NA       NA
##   r[2,1]    0.59   0.61 0.15 0.14   0.30   0.80 1.00      961     1029
##   r[1,2]    0.59   0.61 0.15 0.14   0.30   0.80 1.00      961     1029
##   r[2,2]    1.00   1.00 0.00 0.00   1.00   1.00   NA       NA       NA
##   cv[1,1]   1.19   1.10 0.43 0.36   0.68   1.98 1.00     1360     1029
##   cv[2,1]   0.66   0.60 0.31 0.25   0.26   1.25 1.00      958      816
##   cv[1,2]   0.66   0.60 0.31 0.25   0.26   1.25 1.00      958      816
##   cv[2,2]   1.02   0.95 0.38 0.32   0.58   1.74 1.00     1355     1211

mdl=cmdstan_model('./ex4-44.stan')

k=2
data=list(N=n,K=k,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -64.0672 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       14      -11.9003   4.65303e-05   0.000682734           1           1       26    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##  lp__       -11.90
##  m[1]        -0.09
##  m[2]        -0.01
##  cov[1,1]     0.95
##  cov[2,1]     0.56
##  cov[1,2]     0.56
##  cov[2,2]     0.80
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##  lp__     -13.93 -13.52 1.83 1.58 -17.36 -11.75 1.00      561      792
##  m[1]      -0.09  -0.09 0.28 0.26  -0.54   0.35 1.01      839      806
##  m[2]      -0.02  -0.03 0.24 0.23  -0.41   0.39 1.01      904     1006
##  cov[1,1]   1.46   1.31 0.62 0.49   0.74   2.62 1.00     1262     1064
##  cov[2,1]   0.87   0.75 0.49 0.36   0.32   1.81 1.01     1185      916
##  cov[1,2]   0.87   0.75 0.49 0.36   0.32   1.81 1.01     1185      916
##  cov[2,2]   1.22   1.09 0.53 0.40   0.64   2.21 1.00     1269      990